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ABSTRACT 


A  computer  program  is  described  which  can  directly  determine 
the  composition  which  yields  maximum  impulse  for  multicomponent  propellant 
mixtures . 


The  method  has  been  coded  for  IBM  709  and  7090  computers  and 
has  been  demonstrated  for  systems  containing  up  to  four  components.  The 
mathematics  have  been  determined  so  that  the  technique  is  applicable  to 
systems  containing  up  to  ten  components  but,  thus  far,  it  has  only  been 
applied  to  systems  containing  two  to  five  components. 

The  computation  proceeds  directly  to  the  optimum  point; 
consequently,  an  economy  of  machine  time  over  conventional  procedures 
is  realized.  The  program  can  be  used  in  conjunction  with  any  accurate 
performance  computational  program. 

The  final  report  is  made  up  of  two  volumes:  Volume  I  describes 
the  mathematical  development  and  procedures  adopted  for  carrying  out  the 
optimization  process  and  shows  computer  program  flow  charts;  Volume  II 
presents  the  results  obtained  when  the  program  was  applied  to  some  multi- 
component  systems  of  current  interest  and  discusses  some  interesting 
aspects  of  impulse  surfaces. 
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A  computer  program  has  been  developed  which  can  directly  deter¬ 
mine  rapidly  and  economically  the  propellant  composition  which  yields 
maximum  impulse.  The  method  has  been  coded  for  IBM  709  and  IBM  7090 
computers  and  has  been  demonstrated  for  systems  containing  up  to  four 
components.  The  mathematics  have  been  determined  so  that  the  technique 
is  applicable  to  systems  containing  up  to  ten  components,  but  such  complex 
multicomponent  systems  have  not  been  tested. 

In  the  past,  when  binary  systems  have  been  considered,  the  pro¬ 
cedure  for  finding  the  composition  yielding  maximum  impulse  has  consisted 
of  selecting  weight  ratios  of  the  two  ingredients  such  that  certain  major 
product  species  would  be  formed.  Computations  of  performance  were  made 
of  a  number  of  such  mixture  ratios  and  a  curve  drawn  through  the  points 
to  determine  the  maximum.  Usually,  five  or  more  such  computations,  depend¬ 
ing  upon  the  accuracy  desired,  yielded  the  composition  of  maximum  impulse. 
If  really  exact  determination  of  the  theoretical  maximum  was  desired,  more 
calculations  were  needed,  particularly  if  the  performance  curve  was  irreg¬ 
ular. 


For  three-component  systems,  the  most  convenient  and  reliable 
representation  is  the  use  of  triangular  diagrams.  As  an  alternative  to 
triangular  plots,  Cartesian  coordinates  can  also  be  used,  but  these  can 
lead  to  errors  and  omissions  and  are  not  as  clear  cut.  Using  the  tri¬ 
angular  plot  technique,  a  number  of  compositions  are  again  arbitrarily 
chosen  and  plotted.  When  sufficient  numbers  of  impulse  values  are  de¬ 
termined,  at  least  12  in  the  best  possible  cases  and  usually  many  more 
for  the  average  system,  constant  impulse  contour  lines  are  estimated  and 
drawn  in  the  plot.  Again  the  accuracy  of  the  final  result  is  a  function 
of  the  number  of  computations  made.  A  great  many  more  points  are  usually 
required  for  ternary  systems  than  for  binary  systems,  inasmuch  as  it  is 
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difficult  to  guess  where  the  optimum-performing  composition  will  be.  For 
four  or  more  component  systems,  the  number  of  calculations  required  be¬ 
comes  inordinately  large,  and  representation  of  the  results  and  attain¬ 
ment  of  the  true  maximum  becomes  increasingly  difficult. 

The  computer  program  described  in  this  report  minimizes  the 
number  of  performance  calculations  required  to  determine  the  composition 
of  maximum  impulse.  Valuable  computer  time  is  saved  by  using  the  pre¬ 
viously  calculated  composition,  or  some  close  approximation  to  it,  as  a 
starting  composition  in  each  step  of  the  optimization  process.  In  addi¬ 
tion  to  economizing  machine  time  and  directly  determining  the  optimum 
composition,  the  program  has  been  set  up  to  allow  restriction  of  the  in¬ 
gredients  to  certain  ranges  of  values  or  ratios,  as  desired.  Thus,  for 
a  solid  propellant  system  where  the  presence  of  some  minimum  amount  of 
binder  is  necessary,  even  though  it  degrades  impulse,  the  binder  content 
is  not  allowed  to  drop  below  a  certain  prespecified  value.  Similarly, 
the  ratio  of  oxidizer  to  fuel  can  be  maintained  while  the  binder  content 
can  be  varied,  etc.  These  additions  to  the  program  take  into  account 
practical  considerations  which  cannot  be  ignored  even  in  theoretical  work. 

The  mathematical  details  and  the  programing  of  the  optimization 
procedure  are  described  in  Sections  2  through  7  of  Volume  I.  Some  addi¬ 
tional  programing  to  improve  the  efficiency  of  the  program  is  suggested 
in  Section  8  of  Volume  I.  Volume  II  presents  the  results  obtained  when 
the  program  was  applied  to  some  multicomponent  systems  of  current 
interest . 
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SECTION  2 

MATHEMATICAL  APPROACH  TO  THE  PROBLEM 


2.1  GENERAL  CONSIDERATIONS 

Due  to  the  complexity  of  the  relation  between  the  reactants  in¬ 
volved  in  a  rocket  motor  performance  calculation  and  the  specific  impulse 
produced  by  the  combustion  of  these  reactants,  it  was  felt  that  the  only 
fruitful  approach  to  the  optimization  of  specific  impulse  would  be  one  of 
numerical  nature.  In  the  course  of  developing  the  current  optimization 
program,  two  such  approaches  were  attempted.  These  were: 

(1)  approximation  of  the  impulse  surface  by  a  second 
order  surface,  and 

(2)  a  gradient  approach. 

In  each  of  these  methods  an  initial  composition  of  reactants  is 
chosen,  and  by  considering  certain  properties  of  the  specific  impulse 
function  in  the  neighborhood  of  this  initial  point,  a  new  point  is  deter¬ 
mined.  Hopefully,  the  specific  impulse  at  this  new  point  will  be  greater 
than  that  at  the  old.  Let  us  denote  the  old  and  new  points  by  y  and  z. 
respectively,  and  the  specific  impulse  at  some  arbitrary  point  x  by  I(x). 
Further,  we  will  assume  that  if  y  is  not  the  point  with  maximum  impulse, 
then  the  application  of  our  process  to  y  will  produce  a  z  such  that 
I(z)>I(y).  Clearly,  by  using  z  as  a  new  y  the  process  can  be  repeated 
with  perhaps  a  further  increase  in  I(z).  This  iteration  can  be  continued 
until  no  further  improvement  is  possible.  * 

Note  that  this  procedure  embodies  the  assumption  that  there  is 
only  one  relative  maximum  of  I(x)  in  the  domain  of  x  where  x  is  a  closed 
set.  For  functions  I(x)  that  do  not  satisfy  this  condition  the  procedure 
will  determine  some  point,  say  x* ,  at  which  a  relative  maximum  of  I(x) 
occurs;  however,  I(x*)  will  not  necessarily  be  the  absolute  maximum  of  I(x). 
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2.2  SELECTION  OF  INDEPENDENT  VARIABLES 

The  point  x  must  represent  the  reactants  involved  in  a  par¬ 
ticular  specific  impulse  calculation.  For  instance,  the  i*-*1  coordinate 
of  x  (i.e.,  x^)  could  be  the  mass,  or  perhaps  the  number  of  moles  of  the 
ifck  reactant.  With  either  of  these  definitions  of  x,  however,  we  would 
have  that  I(x)  =  I(kx)  where  k  is  any  positive  constant.  This  follows, 
since  I(x)  is  a  specific  quantity,  i.e.,  independent  of  total  amount,  and 
x  and  kx  represent  mixtures  of  reactants  of  the  same  relative  amounts. 

This  situation,  that  is, where  an  infinity  of  points  represents  a  single 
composition,  has  obvious  computational  disadvantages.  The  basic  problem 
is,  that  for  mixtures  of  n  reactants,  specific  impulse  is  a  function  of 
only  n-1  independent  variables.  For  instance,  we  could  choose  the  ratios 
of  the  amounts  of  the  first  n-1  reactants  to  that  of  the  last  one  as  the 
n-1  coordinates  of  x.  Another  possibility  would  be  to  represent  the 
amounts  of  each  of  the  n  reactants  as  a  coordinate  of  x  but  require  that 

n 

the  total  amount,  that  is  constant.  If  this  constant  is  one, 

i=l  1 

then  the  x^  become  mass  or  mole  fractions. 
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SECTION  3 

OPTIMIZATION  USING  A  QUADRATIC  APPROXIMATION 


The  initial  approach  to  the  optimization  problem  used  a 
quadratic  approximation  to  I(x)  at  y.  This  attempt  was  not  success¬ 
ful,  and  therefore  will  be  treated  only  briefly. 

Let  xi  represent  the  mass  fraction  of  the  ith  reactant  in  the 
mixture.  The  quadratic  approximation  to  I(x)  given  as  a  truncated  Taylor 
Series  about  the  point  y  is 


QM  -  iw+Lje 

i=l  i 


A  r3-  ^  ^2I 

Ai+  f-j  ^x.^x. 
i=l  j=l  a  iCT  j 


x=y 


x=y 


A.  A. 

i  j 


(l) 


where  A.  =  x. -y . . 

i  i  i 

Since  the  x.  are  mass  fractions,  we  have 

i 

x .  =  1 

l 


(2) 


We  will  take  the  new  point  z  to  be  that  point  which  causes  Q(x)  to  be  an 
extremum  subject  to  the  constraint  given  by  equation  (2) .  Applying  the 
Lagrange  Multiplier  technique,  we  are  led  to  the  following  equations: 


A  -  5  1 

Ai  ^  x. 

x=y  j 


j=l,2, 


(3) 


x=y 


where  X  is  a  Lagrange  Multiplier.  Equation  (2)  must  be  satisfied  for  both 
the  old  and  new  points.  It  follows  that 
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A.=0  (4) 

The  partial  derivatives  in  equations  (3)  must,  of  course,  be 
evaluated  numerically.  Equations  (3)  and  (4)  are  a  set  of  n+1  linear 
equations  in  the  unknowns  A  and  the  The  solution  to  this  system 

gives  corrections  which,  when  applied  to  the  point  y,  give  the  new 
point  z.  Thus 

Zi  =  yi+  Ai  i=l»2, • • .n  (5) 

The  algorithm  just  described  makes  the  tacit  assumption  that 
I(x)  can  be  represented  reasonably  well  by  a  second  order  surface  in 
n  dimensional  space.  For  many  systems,  however,  this  is  not  the  case. 
This  can  readily  be  surmised  from  an  examination  of  the  projection  of  the 
impulse  surface  in  various  planes.  This  is  discussed  in  more  detail  in 
the  section  titled  "Characteristics  of  Some  Impulse  Surfaces". 

Due  to  the  poor  approximation  obtained  by  using  second  order 
surfaces,  the  relation  I(z) >I(y)  would  not  in  general  hold.  Therefore, 
this  approacn  to  the  problem  was  abandoned. 
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SECTION  4 

OPTIMIZATION  USING  THE  GRADIENT  APPROACH 


The  method  that  was  actually  ujed  in  the  computer  program  is 
a  modification  of  one  described  by  Curry  .  Other  modifications  of  this 
procedure  have  been  proposed  in  the  literature.  Box's  method  seems  to 
be  an  excellent  one;  however,  it  is  not  easily  mechanized  for  machine 
computation.  The  method  used  by  Booth  for  the  solution  of  linear 
equations  is  clearly  not  adequate  for  the  optimization  of  systems 
whose  specific  impulse  function  is  extremely  rugged.  The  method  des¬ 
cribed  herein  is  easily  mechanized  and  seems  to  be  sufficiently  power¬ 
ful  to  handle  up  to  four  component  systems  with  no  difficulty.  The 
following  discussion  includes  the  option  of  imposing  linear  constraints 
on  the  reactants,  although  this  option  has  not  been  included  in  the  com¬ 
puter  program.  With  the  imposition  of  m  linear  constraints,  the  pro¬ 
cedure  should  easily  handle  systems  of  4+m  components. 

Let  x^  be  the  mass  fraction  of  itb  reactant  of  the  system*. 

As  before,  we  have  the  restriction  that 

n 

3L  *-=i  (6) 

i=l 

Now  suppose  that  we  impose  m  linear  constraints  on  the  x.,,  to  wit: 

aikXi=bk  k=l,2,...m  (7) 

i=l 


*An  optimization  program  was  also  written  using  mass  fraction  ratios  as 
the  components  of  x.  This  procedure  was  somewhat  inferior  to  the  method 
that  was  adopted  for  the  existing  program. 
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Note  that  equation  (6)  is  of  the  same  form  as  equation  (7) .  We  shall 
include  them  altogether  by  writing 

n 

^  a  .  x.=b  k=0,l,2,...m  (8) 

7  i  IK  1  K. 

1=1 


where  bA  and  all  the  a.A  have  the  value  1. 
u  iC) 


Let  us  solve  equation  (8)  for  x^  ,  i»l , 2 . . . m+1  in  terms  of  the 


remaining  x..  Thus 

l 


310  a20 


i  all  a21 


alm  a2m 


‘  nrt-1,0  1 


‘  am+l , 1  |  !  X2 


'  !  • 

i 

a  ,  !  x 

m+l,m  rrH-11 


n 

V.  Wti 

i=nH-2  i 


b  -  V  a  x. 
1  11  i 

i=m+2 


_n_ 

b  -  \  a.  x . 
m  — lm  i 
i=mf  2 


(9) 


Let  A  be  the  determinant  of  coefficients  in  equation  (9)  and  let  A.  be  A 
with  its  ith  column  replaced  by  the  column  vector  on  the  right  side  of  (9). 
Then,  by  Cramer's  Rule 

A. 

x.  =  -j!  i=l , 2  ,  .  .  .n*l  (10) 

1  A 


where  A  is  a  non-zero  constant  and  the  A^  are  functions  of  x^,  i=n*2  ,n*3 .  .  . n 
only.  (The  requirement  that  A^O  is  equivalent  to  requiring  the  constraints 
given  by  equation  (8)  to  be  independent.)  It  follows  that  specific  impulse, 
I  ,  is  a  function  of  only  n-m-1  of  the  reactants;  that  is 


I 

sp 


I(Xmf2’Xn*3’-'-Xn) 


(11) 


Wc  will  denote  the  quantity  n-m-1  as  the  number  of  degrees  of  freedom  of 
the  system. 

In  addition  to  the  constraints  given  by  equation  (8),  we  must  re¬ 
quire  that  the  mass  fraction,  x^,  of  each  reactant  be  non-negative.  It  is 
no  move  difficult  to  require  that  x^  be  not  less  than  an  arbitrary  positive 
constant,  e^,  consistent  with  equation  (8).  For  the  sake  of  generality, 
we  shall  follow  this  course.  Thus,  for  the  n-m-1  independent  variables, 
we  have 
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x.  >  e.  i=mf2  ,mf3.  .  .n  (12) 

and  for  the  remaining  variables  from  equation  (10) 

A.>Aei  i=l,2,...mfl  (13) 

We  summarize  the  results  thus  far  obtained  as  follows:  We  wish 
to  find  the  numbers  x^ ,  i=mf2 ,mf 3, . . .n  subject  to  the  inequalities  (12) 
and  (13)  that  will  cause  ISp,  given  by  the  expression  (11),  to  assume  a 
maximum  value.  (Note  that  once  the  for  i=mf  2  ,irri-3 ,  .  .  .n  are  fixed,  the 

remaining  x^  can  be  obtained  from  equation  (9) . 

The  procedure  at  this  point  as  described  in  Reference  1  would 
be  to  start  at  some  point  y=(ymf2 >ynH-3 > • • • Yn)  and  proceed  in  the  direction 
of  steepest  ascent,  i.e.,  in  the  direction  of  grad  I  ,  to  a  point  z  such 
that  I(z)> I(x)  for  all  points  x  on  the  line  through  y  with  direction 
grad  ISp.  At  this  time,  the  point  z  would  be  treated  as  new  y  and  the 
process  repeated  until  convergence  was  attained.  Such  a  process  is  illus¬ 
trated  in  Figure  1  for  a  system  with  two  degrees  of  freedom.  The  dotted 
line  composed  of  straight  segments  represents  the  path  taken  by  the  iter¬ 
ation.  The  curved  lines  represent  level  curves  of  ISp .  It  is  easily 
shown  that  any  two  adjacent  segments  of  the  path  are  such  that  the  first 
is  tangent  to  a  level  curve  of  ISp  at  the  point  of  intersection  of  the  two 
segments  while  the  second  is  normal  to  the  same  curve  at  the  same  point. 

Figure  1  illustrates  a  surface  whose  maximum  value  is  easily 
found  by  applying  the  algorithm  just  described.  Unfortunately,  the  sit¬ 
uation  is  not  always  this  simple;  for  example,  an  examination  of  Figure  2 
would  indicate  that  very  many  steps  would  be  necessary  in  order  to  locate 
the  peak  for  that  system.  It  would  seem  then  that  we  should  not  be  re¬ 
stricted  to  travel  onlv  in  the  direction  of  grad  I 

sp 

Consider  a  typical  step  in  an  optimization  path  for  a  system 
with  two  degrees  of  freedom,  as  shown  in  Figure  3.  The  points  y'  and  y 
represent  the  previous  and  current  reactant  compositions,  respectively. 

V 1  is  a  unit  vector  in  the  direction  of  the  step  from  y'  to  y,  and  U  is 
a  unit  vector  in  the  direction  of  grad  Igp  at  y.  By  hypothesis,  there  is 
one  and  only  one  relative  maximum  of  ISp  in  dom  x;  suppose  it  occurs  at 
x* .  Denote  the  unit  vector  in  the  direction  from  y  to  x*  by  V,  It  is 
clear  that  if  I  is  continuous  in  dom  x*,  then  U*V>0;  that  is,  U  and  V 


*dom  x  indicates  the  region  in  which  x  is  defined. 
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FI  URE  2.  APPLICATION  OF  THE  GRADIENT  TECHNIQUE  TO  A  LESS  REGULAR 
SYSTEM  THAN  THAT  SHOWN  IN  FIGURE  1. 
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FIGURE  3.  MODIFICATION  OF  THE  GRADIENT  TECHNIQUE 
USING  THE  "T  PROCEDURE" 


S8488 
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will  both  lie  on  the  same  side  of  V'.  It  follows  that  we  can  represent 
V  in  the  form 

V  =  Vl-t2  U+tV'  -l<t<l  (14) 


For  systems  with  more  than  two  degrees  of  freedom,  the  point  x*  will  in 
general  not  be  the  optimum  point  for  the  system,  but  will  be  the  optimum 
point  in  the  plane  defined  by  U  and  V'. 


Suppose  we  allow  the  variable  t  in  equation  (14)  to  take  on  any 
value  between  -1  and  1.  Then  any  point  x,  in  the  half  plane  on  the  side 
of  V'  toward  U,  is  given  by 


x 


=  y+ 


>V  =  y+X  [Vl-t2  U+tV  '1 


(15) 


where 


X 


is  the  distance  of  x  from  y  and  U  has  components 


(16) 


Since  U,  V'  and  y  are  fixed  at  any  point  in  the  iteration,  x  is  a  function 
of  only  /(and  t ..  We  emphasize  this  relation  by  rewriting  equation  (11) 
as  follows: 


i  =  k  (i7) 

bp 

For  a  specified  value  of  t,  ISp  is  a  function  of^only.  Denote  this  function 
as  It(/\).  Let  be  the  value  of  that  maximizes  Ip(^).  Clearly,  ^,v 
will  depend  on  t.  Let  us  restrict  the  ^ in  equation  (17)  to  take  on  only  the 
values  Then 


I  =  K  ^*(t),t)  =  l*(t)  (18) 

that  is,  we  can  consider  ISp  as  a  function  of  t  only.  The  point  at  which 
I*(t)  is  a  maximum  is  clearly  x*. 

The  functions  It(^)  and  I*(t)  can  be  maximized  by  successive 
polynomial  approximations.  This  is  further  discussed  in  the  following 
section . 

The  domain  of  x,  (i.e.,  the  points  representing  acceptable  re¬ 
actant  compositions)  is  defined  by  the  inequalities  (12)  and  (13).  Consider 
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the  equations  formed  from  (12)  and  (13)  by  replacing  the  inequalities 


by  equalities.  Thus, 

x.  =  e . 

i  l 

i=mf 2 ,mf 3 , . . .n 

(19) 

> 

H* 

II 

£ 

h*» 

i=l ,2 , . . .mil 

(20) 

Equations  (19)  and  (20)  represent  planes  in  the  n-m-1  dimen¬ 
sional  x  space.  Dom  x  is  the  region  contained  within  all  the  planes;  that 
is,  if  y  is  some  point  in  dom  x,  then  any  other  point  x  in  dom  x  can  be 
joined  to  y  by  a  line  segment  that  does  not  cross  any  of  the  planes  given 
by  equations  (19)  and  (20) . 

The  above  considerations  make  it  clear  that  Ain  equation  (15) 
may  be  limited  to  certain  values,  say  0<A«£  Amax  where  Amax  *-s  t^e 
value  that  would  cause  x  to  lie  on  the  closest  plane,  given  by  equations 
(19)  and  (20),  in  the  direction  of  V.  Let  us  denote  the  components  of  V 
by  v^ , i=mf2 ,mf 3 , . . .n .  A  line  through  y  in  the  direction  of  V  is  given  by 


x.-y. 
v . 

l 


i=m+2 ,mf 3 , . . .n  (21) 


Now  A  is  the  distance  between  x  and  v.  Denote  the  distance,  along  the 
direction  of  V,  between  y  and  the  jttl  limiting  plane  given  by  equations 
(19)  or  (20)  as  A j •  Then,  from  equation  (19) 


e.-y . 

A  =  J  J 

Aj  v. 

J 


j=m+2 ,m+3 , . . .n 


(22) 


A  .= 
J 


310  a20  •  '  ' 


311  a21  ’ 


(20).  Thus,  for  j=l,2,. 

Vd  -  -  - 

.  .mfl 

3mfl,0 

±2au(V  XjV  ■  ■  ■ 

amfl , 1 

n 

y  a.  (y .+  A  .v.)  .  .  . 
i=m+2  lm  1  J  1 

a 

nrf  i  ,m 

=Ae j  (23) 


-13- 


AERONUTRONIC 


AERONUTRONIC 


B  .= 
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and 


V 


a10  a20  • 


311  a21  ‘ 


3lm  a2m  ‘ 


al0  a20  • 


all  a21  •  * 


alm  a2m 


column  j.  Define 
n 

A  C 

-  .  21  Wi  •  ■  • 

i=nrt-2 

n 

3mf  1 ,0 

-  >  a.  ,y .  .  .  . 

„  ll'i 

3mf  1 , 1 

n 

-  a.  y.  .  .  . 

\  rm  1 
i=mf2 

Q 

m+1  ,m 

a  division  of  ybr<//f(<>/or6e>///Mny, 


u 

i^H-2 

f- 


a..v.  .  .  .  a  ,  . 
lO  1  nH-1,0 


i=nrt-2 


ailVi  ‘  *  *  3m+l , 1 


21«.  v. 

•  To  lm  1 
i=mf  2 


(24) 


(25) 


mf  1  ,m 

where  the  distinct  column  is  column  j .  Then  equation  (23)  can  be  written 


so  that 


-  A-C.=Ae 

j  : 


B . -Ae . 

_J _ 1 

C, 


j=l,2,...mfl  (26) 


Equations  (22)  and  (26)  give  expressions  for  the  distance,  along 
the  direction  of  V,  between  y  and  the  limiting  plane.  The  situation 

is  illustrated  in  Figure  4.  Here  we  have  taken  n=4  and  m=l.  Dorn  x  is  the 
closed  finite  region  bounded  by  the  lines  Ai=Ae^ ,A2=Ae2 1X3=63  and  x^=e^. 

In  this  case  V  points  toward  the  former  two  lines  and  away  from  the  latter 
two.  Thus  here  and^-j,^^<0.  It  is  clear  that  is  the  largest 

step  that  can  be  taken,  in  the  direction  of  V,  starting  from  y,  if  the  point 
x  is  to  remain  in  dom  x.  In  general,  the  largest  step  is  given  by 
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FIGURE  4.  LIMITING  PLANE  RESTRICTIONS  ON  MAXIMUM  STEP  SIZE  IN  DOM  X 
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where  min  pos  stands  for  "the  minimum  positive  value  of"  and 


j-1,2, . . .mfl 


j=m+-2  ,m+3  .  .  .n 
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SECTION  5 

OPTIMIZATION  ON  A  LINE 


The  discussion  thus  far  has  provided  the  basis  for  choosing 
a  region  R,  that  is  a  finite  portion  of  some  half  plane  in  x-space, 
in  which  an  improved  composition  is  to  be  sought.  The  actual  seeking 
out  of  this  improved  point  is  done  by  maximizing  Isp  in  R.  In  order  to 
accomplish  this,  it  is  necessary  to  be  able  to  locate  the  optimum  point 
along  some  line  L  in  R.  The  development  of  an  algorithm  that  will 
efficiently  locate  the  optimum  point  on  L  necessitates  the  assumption 
that  one  and  only  one  relative  maximum  of  Isp  exists  on  L.  This  assump¬ 
tion  does  not  follow  from  our  original  hypothesis;  i.e.,  that  one  and 
only  one  relative  maximum  of  Igp  exists  in  dom  x,  and,  in  fact,  is  not 
generally  valid.  For  example,  see  Figure  5.  Here,  there  is  a  unique 
maximum  in  dom  x;  however,  if  we  restrict  x  to  points  on  L,  then  the 
resulting  function  has  two  relative  maxima*.  None  of  the  cases  that 
have  thus  far  been  run  on  the  computer  have  been  affected  by  the  in¬ 
validity  of  this  assumption.  Since  the  assumption  enables  a  powerful 
tool  to  be  utilized  in  the  optimization  process,  it  will  be  made;  how¬ 
ever,  it  must  be  kept  in  mind  that  this  procedure  could  provide  a  possible 
source  of  difficulty  in  the  optimization  of  a  given  system.  As  will  be 
noted  below,  this  assumption  can,  in  effect,  be  removed  from  the  optimi¬ 
zation  process  by  the  suitable  choice  of  a  program  constant.  This,  of 
course,  should  not  be  done  unless  trouble  arises  in  using  the  recommended 
procedure . 


The  difficulty  arises  in  the  maximization  of  the  function  I*(t)  given  in 
equation  (18).  The  numerical  procedure  used  in  computing  I*(t)  is  such 
that  it  is  possible  for  either  of  the  maxima  indicated  in  Figure  5  to  be 
specified  as  the  optimum  point  on  L.  The  resulting  function  I*(t)  could 
therefore  be  discontinuous. 
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FIGURE  b.  TECHNIQUE  FOR  LOCATING  MAXIMUM  LYING  WITHIN  DOM  X 
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Consider  the  function  defined  by  equation  (17)  for  a  spec¬ 
ified  value  of  t.  The  region  of  definition  for  this  function  is  L, 
the  line  segment  composed  of  the  points  x  given  by  equation  (15)  for 
0  <  X<Amax  where  X  max  *-s  8*-ven  by  equation  (27).  We  have  already 
referred  to  this  function  as  It(A).  It  will  be  assumed  that  one  and 
only  one  relative  maximum  of  It(  A)  exists.  Under  this  assumption, 
there  are  three  possible  situations: 

(a)  It(0)  is  the  relative  maximum, 

(b)  I  (A  )  is  the  relative  maximum,  or 

t  'max 

(c)  a  relative  maximum  exists  for  0  <XcA  . 

max 

The  process  of  locating  the  optimum  A  offers  no  problems  in  the  first 
two  cases.  Accordingly,  we  will  examine  the  third  case. 

Consider  a  graphical  representation  of  It(A)  in  which  the 
ordinate  and  abscissa  of  a  point  represents  It(  A-^)  andAi>  respect¬ 
ively.  Clearly,  we  can  determine  three  values  of  %  ,  say  Ai<  A2^  A3> 
such  that  I.  (  A2)  ^  It(  Ai)  >  i=l,3.  Through  the  corresponding  points 
P^,  P2  and  P3,  let  us  construct  a  parabola,  and  then  determine  the  value 
of  A  corresponding  to  the  vertex  of  the  parabola,  sayAp.  Next,  we  can 
compute  It(  Ap)and  thereby  obtain  Pp=(  Ap » It ( A p) ) •  From  the  four  points 
Pi,  P2.  P3  and  P  we  can  choose  that  with  the  largest  ordinate.  (This, 
of  course,  will  be  either  ?2  or  Pp  • )  Using  this  as  a  new  P£ ,  together 
with  the  adjacent  points  on  either  side  of  it  as  new  P^  and  P3,  we  ob¬ 
tain  a  new  set  of  three  points,  and  the  process  can  be  repeated.  This 
is  illustrated  in  Figure  6.  Here,  the  next  set  P^ ,  ?2>  P3  would  be  the 
old  points  P^ ,  Pp,  P2 . 

The  program  actually  employs  a  process  that  is  a  slight  mod¬ 
ification  of  the  procedure  just  described.  The  trouble  with  the  above 
procedure  is  that  it  can  be  rather  slow  in  converging.  For  example, 
consider  Figure  7.  The  initial  three  points  are  A,  B  and  C.  In  this 
case  the  curve  is  such  that  a  poor  approximation  to  the  optimum  A  is 
obtained  at  each  step.  The  process  yields,  successively,  the  points 
D,  E,  F,  G  and  H,  and,  in  turn,  discards,  successively,  the  points 
B,  C,  D,  E  and  F.  Here,  of  course,  the  troublesome  point  is  A,  and  its 
retention  prevents  good  approximation  to  the  optimum 


In  order  to  alleviate  this  difficulty,  the  current  procedure 
sometimes  determines,  in  addition  to  Pp ,  the  point  Pm=( Am,I ( Am))  where 
Am=^(  A  A3)  •  This  is  done  only  if  Am  considerably  different  than 
both  ^2  and  ^  •  The  criterion  being  used  is  that  if  both 
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then  Pm  is  determined.  The  actual  process 


1. 

Start  with  a  set  of  points 

(P1> 

P  Pi 

2’  3' 

• 

2. 

Compute  and  . 

3. 

0  -  1 

If  .  m-  <0.2,  go  to  8. 

a3~  Ai 

4. 

If  <0.2,  go  to  8. 

A3'  Ai 

5. 

Otherwise,  compute  I^C^^), 

thereby  obtaining  P  . 

6. 

Get  new  set  from 

old 

points 

<PX-P2-P3>Pm> 

as  described  above. 

7. 

Compute  from  new  set. 

8. 

Compute  It(/lp)  thereby  obtaining  P^ . 

9. 

Get  new  set  (P^.P^.P^)  from 

old 

points 

<P1'P2'P3'V 

as  described  above. 

10. 

Go  to  2 . 

Whenever  a  new  set  of  points  is  obtained,  they  are  tested  in  order  to 
determine  if  convergence  has  been  achieved  and  therefore  if  the  iteration 
should  be  terminated. 

The  foregoing  discussion  gives  the  method  used  in  evaluating 
I*(t)  defined  in  equation  (18).  Once  t  is  specified,  the  optimum  A  > 
i.e.,A*(t)>  can  t>e  determined  by  the  above  procedure.  Thus  I*(t)  is 
obtained.  In  finding  the  value  of  t  that  will  maximize  I*(t) ,  a  similar 
procedure  is  used.  In  this  case  the  ten  steps  outlined  above  are  repeated 
except  that  all  It(A)  ate  replaced  with  I ( t )  and  all  ^  are  replaced  with  t. 
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The  determination  of  the  maximum  value  of  I*(t)  is  a  rather 
lengthy  process  when  compared  to  the  computation  involved  in  the  max¬ 
imization  of  It( A).  In  order  to  avoid  unnecessary  computation  and 
still  retain  the  use  of  the  I*(t)  maximization  when  applicable,  the 
following  procedure  is  effected. 

(a)  Grad  Isp  is  determined,  and  IQ( A)  Is  maximized. 

This  computation  provides  ^*(0)  as  defined  in 
conjunction  with  equation  (18) . 

(b)  If  A*(0)^e^  where  e^  is  a  program  constant 
(presently,  we  are  using  0.02),  then  go  to  (a). 

(c)  Otherwise  maximize  I*(t). 

(d)  Go  to  (a) . 

Note  that  if  e  is  taken  to  be  zero,  then  the  maximization  of  I*(t) 
is  completely  bypassed,  and  the  process  reduces  to  the  usual  gradient 
method  as  described  by  Curry,  in  Reference  1.  In  any  case,  where  the 
assumption  of  a  unique  maximum  of  I  on  a  line  causes  computational 
difficulties,  this  procedure  can  bes^followed . 
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SECTION  6 

DETERMINATION  OF  THE  GRADIENT 


(X  X  , .  ..XL  )  be  an  arbitrary  point  in  the  n-m-1 
m-2,  nri-3,  n 


Let  P 

dimensional  x-space.  The  impulse  function,  I 


surface,  I(P),  in  the  n-m  dimensional  space. 


,  can  be  considered  as  a 

?f  P  =  (Y  Y  ...Y  ) 
o  m+2,  mf3,  n 
is  some  specified  point  in  the  x-space,  then  the  tangient  hyperplane 

to  the  impulse  surface  at  P^  is  given  by 


I(P)  -  I(PQ)  - 


II 

z 

i=m+2 


d  I 

'dx. 


(X.  -  Y.) 
1  l 


=  0 


(28) 


The  derivatives  in  equation  (28)  are  the  components  of  grad  ISp  at  p 
Since  they  are  not  directly  obtainable  we  must  use  equation  (28)  in 
order  to  evaluate  them  numerically.  Thus  we  can  choose  points  P.,  j  = 

mf2,  m+3,...n  in  the  neighborhood  of  P  and,  after  evaluating  the  I(P.), 


equation  (28)  yields  n-m-1  linear  equations  in  the  unknowns  — ^-A 

a  X . 


Thus 


I(P.)  -  I(PQ) 


II 

-  y  ^i_ 

i^2  dXi 


(Xi  -  Yi> 


itH-2,  . 


(29) 


where  X^  is  the  ith  coordinate  of  P . .  The  solution  to  equations  (29) 

affords  an  approximation  to  grad  I  at  P  . 

sp  o 


-23- 


AERONUTRONIC 


AERONUTRONIC 


a  division  of  ybrdj%o/ortfo»i/ia>i 


A  simpler  procedure  would  be  to  choose  the  P.  such  that  only 

the  coordinate  of  P.  differs  from  that  of  P  .  In  ilhat  case  the  system 

J  J  o  1 

of  equations  (29)  would  reduce  to 


I(P  )  -  I(P  ) 

— - - - -  i  =  nH-2>m_3,  .  .  .n  (30) 


If  there  is  some  a  priori  knowledge  of  the  nature  of  the  impulse  function 
it  may  be  feasible  to  choose  the  Pj  in  such  a  way  as  to  minimize  the  error 
inherent  in  the  approximate  equations  (29).  However,  in  general  we  do 
not  have  such  knowledge.  Therefore,  we  shall  prefer  the  simpler  equation  (30). 

Let  us  define  the  point  AP^  as  follows 


Al_ 

ax. 


Ap. 

i 


(0,0,0,...AX.,0,0,...0) 


(31) 


where  AXj=  X^  -  is  the  i  coordinate  of  the  right  side  of  equation  (31). 
Then  equation  (30)  becomes 


ai 


axi 


I(Po  +  AP.)  -  I(P0) 
ax! 

i 


m+-2 ,  .  . 


(32) 


Equations  (32)  are  linear  approximations  to  the  derivatives  at  PQ. 

It  has  been  found  that  by  improving  the  accuracy  of  the  derivatives 
at  each  step  the  iteration  to  optimum  I  is  accelerated.  A  more  accurate 
set  of  derivatives  can  usually  be  obtainld  if,  instead  of  using  equations 
(32),  a  quadratic  approximation  is  applied.  If  the  points  chosen  for  the 


quadratic  approximation  to 


approximating  form  is 


AI_ 

c)X 


are  P  -  A  P  , 
o  1 


P  +  Ap,  then  the 
o  1 
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If  equations  (33)  are  used  instead  of  equations  (31),  an  additional  impulse 
must  be  obtained  for  each  i.  The  additional  computation  time,  however,  is 
small  since  AX^  is  chosen  to  be  a  small  number*  and  so  a  good  guess  for  the 

composition  at  -  AP .  is  available.  (The  guess,  of  course,  is  the  com¬ 
position  at  P  . )  The  program  accordingly  uses  equations  (32)  to  obtain 
the  gradient. 


Currently  we  are  using  AX^  =  max  (0.0015625  X^,  0.000625) 
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SECTION  7 

PERFORMANCE  PROGRAM  REQUIREMENTS 


The  derivatives 


^  I 

^Xi 


specified  in  equation  (16)  must  be 


puted  numerically.  In  order  to  determine  these  fairly  accurately, 
small  increments  in  x^  must  be  utilized.  This  is  especially  true 
systems  with  narrow  ridges  in  the  neighborhood  of  such  ridges  and 
system  in  the  vicinity  of  its  peak.  This  requirement  of  small  xi 
ments  necessitates  a  very  precise  specific  impulse  calculation. 


corn- 

rat  her 
for 

for  any 
incre- 


Another  facet  of  the  optimization  procedure  that  requires  a 
precise  evaluation  of  Isp  is  the  one  dimensional  maximization  of  It(A) 
and  also  that  of  I*(t).  In  the  vicinity  of  a  one  dimensional  relative 
maximum  it  is  likely  that  the  approximating  curve  will  be  computed  from 
points  at  close  proximity.  Thus  a  small  deviation  in  the  value  of  Isp 
at  a  point  might  radically  change  the  nature  of  the  approximating  curve 


In  computing  numerical  derivatives  the  computer  program  uses 
increments  in  x^  of  between  0.000625  and  0.0015.  The  program  determines 
I*(t)  ,  i.e.,  the  maximum  value  of  I(.(  A),  to  an  accuracy  of  0.001  in  X 
and  0.001  in  Isp-  These  factors  make  it  imperative  that  the  performance 
program  be  precise  (i.e.,  continuous)  to  within  0.0005  seconds  in  I 

sp 

The  convergence  criteria  mentioned  above  could,  of  course,  be 
somewhat  relaxed.  It  has  been  noted,  however,  that  in  the  rather  prevalent 
case  of  a  surface  containing  a  narrow  ridge,  the  determination  of  the  peak 
impulse  is  accelerated  by  requiring  rather  precise  convergence  to  be  at¬ 
tained  at  intermediate  points  in  the  iteration.  This  is  further  discussed 
in  the  next  section. 
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SECTION  8 

SUGGESTED  ADDITIONAL  WORK 


In  order  to  make  the  program  more  useful  and  also  more 
efficient,  the  following  tasks  should  be  undertaken: 

(1)  The  option  to  impose  linear  constraints  on 
the  system  should  be  programmed.  The  math¬ 
ematical  formulation  of  the  inclusion  of  such 
constraints  has  been  carried  out  as  previously 
indicated . 

(2)  An  improved  method  for  finding  the  peak  im¬ 
pulse  along  a  line  in  mass  fraction  space 
should  be  sought.  It  would  seem  that  a  com¬ 
bination  of  cubic  and  parabolic  approximations 
to  the  peak  impulse  would  be  preferable  to  the 
parabolic  and  mid-point  approximation  pre¬ 
sently  being  used. 

(3)  If  it  is  desired  to  optimize  systems  with  four 
or  more  degrees  of  freedom,  some  consideration 
should  be  given  to  a  generalization  of  the  "t" 
iteration  to  higher  dimensions. 

(4)  Special  tests  for  convergence  should  be  devised. 
In  establishing  convergence  to  the  optimum  pro¬ 
pellant  composition  the  current  computer  programs 
sometimes  take  an  undue  amount  of  time. 
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SECTION  9 
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APPENDIX  D 

SUBROUTINE  FOR  COMPUTING  OPTIMUM  VALUE  OF  A  FUNCTION  OF  ONE  VARIABLE 

The  following  symbols  are  used  in  extrme  charts  1 . through  4: 

XR  -  index  register 

FUNCT  -  The  address  of  the  starting  location  of  the  subroutine 

that  computes  the  function  to  be  extremized.  NOTE:  This 
subroutine  must  preserve  index  Register  1. 

N  -  The  maximum  number  of  times  that  the  function  is  to  be  computed. 

TABLE  -  A  block  of  storage  in  which  all  pertinent  quantities  are  saved 
in  the  extremization  process. 

T  -  The  level  of  the  extremization  (the  optimization  of  I  (  X  )  is 

the  0-th  level  and  the  optimization  of  I*(t)  is  the  1st  level.) 

X  -  The  starting  location  of  a  block  of  three  words  in  which  the 

abscissas  of  the  bracketing  points  are  to  be  found. 

Y  -  Same  as  X  except  the  ordinates  are  stored  here. 

ERRX  -  Address  of  cell  containing  tolerance  for  the  abscissa 

ERRY  -  Address  of  cell  containing  tolerance  for  the  ordinate 

P(x,y)  -  The  point  P  having  abscissa  and  ordinate  x  and  y 

m(subscript)  -  "At  midpoint" 

p(subscript)  -  "At  parabolic  approximation" 

1,  2,  and  3  (subscripts)  -  Leftmost,  middle  and  rightmost  points  in 
the  set  of  three  bracketing  points 
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APPENDIX  D  (Continued) 

EXTRME  CHART  1 


SAVE  XR‘S  2  8  4  IN 
TEMPORARY  LOCATIONS. 
TURN  OFF  OVER  FLOW 
LIGHT. 

PICK  UP  4th  « 
CALLING  SEQ 
XR2.  STORE 
TEMPORARY  L 

IORD  IN 
FUNCT— ► 

N  IN 

DC  AT  ION 

COMPUTE  POSITION  IN 
TABLE  OF  FIRST  LOCATION 
IN  BLOCK  OF  CELLS 
UTILIZED  FOR  THE  Tth 
OPTIMIZATION  LEVEL.  THIS 
QUANTITY  GOES  INTO 
XRI _ _ 


STORE  THE  FOLLOWING 
IN  THE  APPROPRIATE 
CELL  OF  Tth  BLOCK: 
i:  TRA  FUNCT  INSTRUCTION 

2.  N 

3.  XR  S  2  8  4 

I 


PICK  UP  2nd 

OF  CALLING 

P  Z  E  X,  0 

WORD 

SEQUENCE 

.  * 

STORE  COORDINATES  OF 

THE  THREE  POINTS 
BRACKETING  EXTREME 

VALUE  IN  APPROPRIATE 
CELLS  OF  T  th  BLOCK 

/ARE  THE  THREE 
/POINTS  PROPERLY' 
(ORDERED  B  DO 
(THEY  BRACKET 

\the  extreme 

N.  VALUE 


NO 


Iyes 

[  PICK  UP  3rd  WORD 

of! 

CALLING  SEQUENCE  j 

PZE  ERRX  ,  0,  ER  RY 

8 

STORE  TOLERANCES 

IN 

APPROPRIATE  CELLS 
Tth  BLOCK 

OF 

/are  any  V 

(^TOLERANCES  .TO,y- 


YES 


COMPUTE  Xm 


/ 


X 


COMPUTE 

XP® 


IS  EITHER  Xm 
OR  Xp  CLOSE 
TO  X2 


NO 

I-  / 

Xm  TO 

ACCUMULATOR 

h_ 

— v 

0 


NO 


0 - » 

RESTORE 

XR’S  2  8  4 

_ w  error 

”  EX  IT 

_ 

PUT  X  AND  Y 
COORDS  OF  EXTREME 
POINT  IN  MQ  8 
ACCUMULATOR 

EXIT 


1 


7 
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APPENDIX  D  (Continued) 
EXTRME  CHART  2 


WHERE 


ASY|  (X2-X3j  +Y2  (X3  -  X,)  +  Y3  ^X|  -  X2^ 

B  =  Y,  (x*-  A)  +Y2  04"  X?)  +  Y3  (X^-  X*) 
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APPENDIX  D  (Continued) 
EXTRME  CHART  3 

(D  COMPUTE  ORDINATE 
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APPENDIX  D  (Continued) 

EXTRME  CHART  4 

(3)  GET  NEW  SET  OF  POINTS 
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APPENDIX  E  (Continued) 
ADDITIONAL  SUBROUTINES 
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APPENDIX  F 

OPERATING  INSTRUCTIONS  FOR  USE  OF  OPTIMIZATION  SUBROUTINE 


The  propellant  optimization  routine  is  an  open  subroutine. 

Several  parameters  and  a  function  subroutine  must  be  supplied  by  the 
user.  The  parameters  and  their  symbolic  names  are: 

(1)  The  number  of  propellant  ingredients;  KPROP 

(2)  The  mass  fractions  (or  weight  percents)  of  the 
ingredients,  their  chemical  formula  and  heat  of 
formation  at  298°K;  block  starting  at  PROPMF 
(up  to  10  cells) 

(3)  The  gaseous  product  species  mass  fractions;  block 
starting  at  XG  (up  to  100  cells) 

(4)  The  condensed  product  species  mass  fractions;  block 
starting  at  XC  (up  to  20  cells) 

The  function  subroutine  to  be  supplied  is  assumed  to  have  the 
symbolic  entry  IMPULS  (index  register  4  is  set  upon  entry)  and  should 
return  to  the  optimization  routine  by  TRA  1,  4  with  the  function  value 
in  the  accumulator. 

The  specification  of  the  reactant  mass  fractions  allows  the 
user  to  specify  the  point  on  the  impulse  surface  from  which  the  opti¬ 
mization  routine  begins.  The  product  species  mass  fractions  are  used  to 
permit  the  optimization  routine  to  provide  the  IMPULS  subroutine  with 
good  estimates  in  order  to  minimize  the  computation  time  required  to 
calculate  specific  impulse.  The  use  of  the  latter  capability  is  not 
required. 

With  the  completion  of  initialization,  control  is  transferred 
to  the  optimization  routine  with  TSX  OPTMUM,  4.  Specific  impulse  is 
calculated  for  the  original  fuel  composition  and  that  is  regarded  as  the 
first  base  point.  The  derivatives  of  impulse  with  respect  to  changes  in 
fuel  composition  are  determined  numerically.  2n-2  impulse  calculations 
(n  =  the  number  of  reactant  species)  are  required  for  derivative  calculations 
at  a  base  point.  These  derivatives  are  used  to  define  the  line  along  which  a 
maximum  impulse  is  sought.  If  the  "t"  iteration  is  used  (only  after  at  least 
two  base  points  are  established)  the  line  is  defined  as  a  linear  combination  of 
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the  lines  established  by  the  numerical  derivatives  at  two  successive 
base  points.  After  the  maximum  impulse  along  the  line  (or  among  the 
maxima  of  several  lines  if  the  t  iteration  is  used)  is  found,  the  point 
at  which  it  occurs  is  used  as  the  next  base  point  and  the  computation 
proceeds  from  there  as  it  did  from  the  previous  base  point.  Conver¬ 
gence  is  identified  in  one  of  two  ways: 

(1)  The  impulse  values  at  three  successive  base  points 
are  within  an  interval  of  one  second, 

(2)  The  use  of  the  "t"  iteration  offers  no  improvement 
to  the  value  of  specific  impulse  at  the  current 
base  point. 

If  in  the  search  for  maximum  specific  impulse  a  reactant  species 
concentration  tends  to  zero,  an  attempt  is  made  to  allow  the  iteration  to 
proceed  with  only  a  small  amount  of  that  species  present.  If  the  attempted 
elimination  of  the  species  is  persistent  the  computation  is  terminated 

After  convergence  or  termination  of  the  optimization, control 
is  transferred  to  the  instruction  following  the  transfer  to  the  optimiza¬ 
tion  subroutine. 


